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Abstract. 

By solving the three-dimensional Gross-Pitaevskii equation we generate two- 
dimensional axially-symmetric and anisotropic dipolar Bose-Einstein condensate bright 
solitons, for repulsive atomic interaction, stabilized by only a weak one-dimensional 
optical lattice (OL) aligned along and perpendicular, respectively, to the dipole 
polarization direction. In the former case vortex solitons can also be created. We 
show that it is possible to make a stable array of small interacting axially-symmetric 
dipolar solitons put on alternate OL sites. Further, we demonstrate the elastic nature 
of the collision of two such solitons. 
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A bright soliton is a self-reinforcing solitary wave that maintains its shape, while 
traveling at constant speed, due to a cancellation of nonlinear attraction and dispersive 
effects. Integrable solitons without any external trap or intervention for cubic 
nonlinearity exist only in one dimension (ID). Experimentally, bright matter- wave 
solitons and soliton trains were created in a quasi-lD Bose-Einstein condensate (BEC) of 
^Li [H E] and ^^Rb atoms [3j by turning the atomic interaction attractive from repulsive 
using a Feshbach resonance (FR) [4J and employing a transverse trap. 

Although, the normal three-dimensional (3D) BEC soliton [5j is of great interest, 
such a BEC has only short-range attraction which makes it vulnerable against collapse. 
Physical systems are stable due to a peculiar nature of interaction among its constituents 
(atoms, molecules and nuclei), e.g., short-range repulsion and long-range attraction. 
Lately, BEC of ^^Cr O [7j and ^^^Dy [HI [9] atoms with a large long-range dipolar 
interaction has been observed. Also, experimental tuning of the long-range dipolar 
interaction by means of rapidly rotating orienting fields [lOj as well as of the short- 
range atomic interaction using a FR [4J are completely under control. This engineering 
of the atomic and dipolar interactions makes the dipolar BEC (DBEC) an interesting 
system for the formation of soliton [HI [121 US]. The long-range anisotropic dipolar 
interaction is attractive in some directions and repulsive in others. If it were attractive 
in all directions, stable robust 3D DBEC solitons, corresponding to a minimum in energy 
functional, would naturally be formed for repulsive short-range atomic and attractive 
long-range dipolar interactions [14j . 

Normal dipolar interaction leads to attraction along the polarization z direction and 
repulsion along transverse directions. It is possible to have the opposite by tuning the 
dipole interaction to "negative" values by orienting fields and this set-up was used 
in some studies [T3l ITS] . For normal dipolar interaction, an anisotropic two-dimensional 
(2D) soliton can be obtained for repulsive short-range atomic interaction if a weak 
OL is placed along y axis, perpendicular to the polarization direction z, to overcome 
the dipolar repulsion in transverse directions. For dipolar interaction tuned to negative 
values [To] , axially-symmetric 2D bright and vortex solitons can be obtained for repulsive 
short-range atomic interaction if a weak OL is placed along z axis to overcome the 
dipolar repulsion in that direction in this sign-changed setting. In all cases, the dipolar 
repulsion is weak and we do not need any trap in other directions to stabilize a soliton. 
Such 2D solitons cannot be stabilized without the dipolar interaction [T3] . 

We present a linear stability analysis for the axially-symmetric soliton [16j. We 
study the 2D DBEC solitons using the numerical and Lagrangian variational analysis 
of the 3D Gross-Pitaevskii (CP) equation. The effective Lagrangian of the variational 
analysis has the same structure as that of a generalized classical dynamical system with 
two degrees of freedom. We find that stable (referred to as "center" as it corresponds to 
a stable periodic orbit around a center in phase space) and unstable stationary (called 
"saddle" as it corresponds to a saddle point in energy) states appear and disappear 
through the mechanism of saddle-center bifurcation p7]. 

There have been studies of 2D DBEC solitons with strong harmonic traps along 
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y [TT] or z [13] axis and of ID DBEC solitons under transverse harmonic trap [T2]. 
The present sohtons confined by only a weak OL along y or z axis, respectively, are 
distinct. The previous studies [Ill[13l[l5] will essentially have an approximate Gaussian 
density distribution along the infinite trap direction, whereas the present solitons will 
have an exponential density distribution due to weak finite traps in these directions. 
More interestingly, an OL simulates the periodic electron-atom potential in a solid and 
the study of solitons in an OL is also of interest in condensed-matter physics [181 ttS] • 
We show that a new type of stable interacting ID array of solitons can be formed in 
3D space when tiny axially-symmetric interacting DBEC bright solitons are placed on 
alternate sites of the OL. However, if the solitons are placed on all sites of the OL, the 
array is destroyed due to strong long-range dipolar interaction among its constituents. 
Statics and dynamics of such periodic array of tiny droplets of dipolar matter are of 
concern in condensed matter physics [19], as they simulate many problems of general 
interest, such as, a periodic linear array of tiny magnets. Polarized droplets of ^^Cr and 
^^^Dy have permanent magnetic dipole moment. 

In a repulsive BEG on 3D OL, gap solitons having negative effective mass 
responsible for attraction, with the chemical potential lying in the band-gap, can be 
made [201 [2T] . The present solitons on ID OL, free to move in the transverse plane are 
bright, and not gap, solitons. 

We consider a DBEG of atoms, each of mass m, using the GP equation: [6] 

~ + Vl)l + 9\<P\' + F cP{r,t) (1) 

with g = 47raN, F = / [/,,(r - r')|(/)(r', r = {x,y,z} = {p,z}, Vg^ = 

— T4cos(2z) or —VyC08{2y) is the weak OL for stabilizing the soliton, Udd{^) = 
9dd{^ — ^cos'^ S)/R^^gdd = Scidd^^^^ R = F — r', normalization / (/)(r)^rfr = 1, a the 
scattering length, 9 the angle between R and z, add = fiofi^m/{127rh^) the strength of 
dipolar interaction, the (magnetic) dipole moment of an atom, and /xq the permeability 
of free space. The parameter a {1 > a > —1/2) can be tuned by a rapidly rotating 
magnetic field allowing the change of the sign of dipole interaction. In ([T]), length is 
measured in units of Iq = A/(27r), time t in units of to = mll/h^ Vy^Vz^ and energy in 
units of 2Eji^ where Eji = /{2m\^) is recoil energy, with A the OL wave length. 

First we consider the axially-symmetric soliton for Vq^ = — cos(2z) and a < 0. 
In this case the Lagrangian density of Q is [221 ES] 

C=\i (00* - 4>*4>t) + ^|V0|2 + 27raiV|0r + V(^lI'^I' 

+ \N\<t>? juM{r-v')W)\Hv'. (2) 

For a variational study we use the Gaussian ansatz [221123!: 0(i", ^) = exp(— p^/2w^ — 
/2w1 +i'-fp^ + if3z'^)/{wp^/w^7^^^'^) where Wp and are time-dependent widths and 7 
and /3 are time-dependent chirps. The effective Lagrangian L (per particle) is 

L= I Cdv = 



' dt 
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Figure 1. (Color online) (a) Numerical 3D contour of an axially- symmetric soliton 
for ^ = 50 and gdd = —15 on OL Vq£ = — 2cos(2z). (b) The same for an anisotropic 
soliton with ^ = 50 and gdd = 20 on OL Vq^ = —2cos{2y) and that for a vortex 
soliton with g = 5 and gdd = —9 on OL Vq^ = — 2cos(2z) at t = (c) and (d) 70. 
The density |^(r)p on the contour is 0.001. 



+ ^kin + ^trap + ^int, (3) 

with kinetic, trap, and interaction energies given, respectively, by E'km = {l/2w'^ + 
l/4:wl),Etr^^ = -V^exp{-wl), Ei^t = N[a - addf{i^)]/{V^wlw^), where /(/^) = 
[1 + 2/^2 _ 3^2 d{K)]/{l - K^), d{K) = (atanhyr^^)/yr^^, k = Wp/w,. The Euler- 
Lagrange equations for parameters can be used to obtain the following 

equations of the widths for the dynamics of the DBEC state 

Wp = ^ + -7=^ — 2a - adde[i^)\ , (4) 



t^p V 27r w^^w^ 

Wz = + -y=^^ [a - addh[n)\ -— -, (5) 

withe(/^) = [2-7H.^-4.^^+9^^d{H)]/{l-H.^f,h{H) = [l + WK^-2K'^-9K^d{K)]/{l-K^y. 
The widths of a stationary soliton of energy E = E'kin + £'trap + ^int are obtained by 
solving ^ and ^ for Wp = Wz = 0. 

To obtain a quantized vortex of unit angular momentum h] around z axis, we 
introduce a phase (equal to the azimuthal angle) in wave function [2l] . This procedure 
introduces a centrifugal term l/[2(x^ + y'^)] in the GP equation for a vortex and we 
adopt this method to study an axially-symmetric vortex soliton on a ID OL along z 
axis for the dipolar interaction tuned to negative values. 

For the anisotropic 2D soliton on OL Vq^ = —Vy cos(2y) with a > 0, we consider a 
minimization of energy E for a soliton using the Gaussian ansatz (j){r) = exp{—x^ /2wl — 
y^/2wl - z^/2wl)/Qw,WyW,iT^'^), with E^^^ = l/4wl + 1/4^2 ^ 1/4^2^ ^^^^^ = 
-y^exp(-^^), Eir,t = [a + adds{K, ky) - add\l\\f2n\w^WyW^ and 

s[K,ky)= 6) 

'l + {^l-l)u\/l + {^l-l)u^ 



where = uJx/uJz, i^y = ^jy/^Jz- 

We perform numerical simulation of the 3D GP equation ([T]) using the split-step 
Crank-Nicolson method [25^. The dipolar term is treated by fast Fourier transformation 
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Table 1. Numerical (n) and variational {v) energy and rms sizes E , (x) , (y) ^ (z) of 
bright (br) and vortex (vor) solitons. Experimental parameters add (15ao for ^^Cr 
and 130ao for ^^^Dy), a, N for realizing these solitons are given for scattering length 
a = 5ao and wavelength A = 10000 A. 
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Q-dd (co) 


a 
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E 


{x) 


{y) 




n, br 


50 


-15 


15 


-0.418 


2400 


-0.836 


2.00 


2.00 


0.450 


V, br 


50 


-15 


130 


-0.0482 


2400 


-0.814 


2.086 


2.086 


0.423 


n, br 


50 


20 


15 


0.557 


2400 


-0.752 


0.533 


1.38 


3.65 


n, vor 


5 


-9 


130 


-0.289 


240 


-0.744 


4.02 


4.02 


0.663 



[22] . The error of the reported numerical results is less than 1 %. We present in figure 
[l] (a) the 3D contour of the axially-symmetric bright soliton for g = 50, gdd = — 15 and 
Vq^ = — 2cos(2z). In figure [l] (b), we show the anisotropic bright soliton for g = 50, 
gdd = 20 and Vq^ = — 2cos(2y). For the anisotropic soliton, the numerical energy is 
—0.752 in agreement with the energy —0.739 obtained from the minimization in ([6]). 
The anisotropy in the x — z plane in figure [l] (b) is due to dipolar interaction. In figure [l] 
(c) we show an axially-symmetric vortex soliton on ID OL, Vq^ = —2 cos(2z) for g = 5 
and gdd = — 9. A relatively large \gdd\ is needed to overcome the centrifugal barrier 
and stabilize a vortex soliton. The bright solitons of figures |l| (a) and (b) are stable 
in real-time propagation. However, the vortex soliton with the parameters of figure [T] 
(c) suffers from transverse instability at large times {t > 70), which eventually leads to 
its destruction [26]. The snapshot of the vortex soliton after real-time propagation at 
t = 70 in figure [l] (d) does not, however, show any distortion or sign of instability. The 
numerical energy and root-mean-square (rms) sizes of the solitons of figures [l] are shown 
in Table I with variational results in the axially-symmetric case. In this table we also 
show the parameters add^ and for these solitons for ^^Cr and ^^"^Dy stabilized by a 
laser of wavelength A = 10, 000 A, and atomic scattering length a = 5ao obtained using 
a Feshbach resonance. 

For the axially-symmetric bright soliton, we have a conservative system with two 
degrees of freedom Wp and with Lagrangian The stable state appears and 

disappears by saddle-center bifurcation [17] as is increased as shown in figure [2] (a) 
for g = 50 and gdd = —15, where the unstable stationary (Mi and M2) and stable (S) 
states are shown in the Wp versus plot. For small Vz{< 0.4278) there exists only the 
unstable stationary state Mi. At Vz = 0.4278 and Wp ^ 10 a stable (S) and a unstable 
stationary (M2) state appear "out of nothing" by saddle-center bifurcation. With further 
increase of V^, the center S comes towards the saddle Mi and the two disappear "to 
nothing" by a reverse (sub-critical) saddle-center bifurcation at = 12.2, whereas the 
state M2 moves towards infinity [17j. We show the equal-energy variational contours in 
the Wz versus Wp phase plot for different Vz in figures [2] (b) — (h) , where the positions of 
the stable and unstable stationary states are also shown. Figures [2] (d) — (g) show close 
ups of the appearance and disappearance of the state S by saddle-center bifurcations. 
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Figure 2. (Color online) (a) Bifurcation diagram showing saddle-center bifurcation 
involving the stable (center, S) and unstable stationary states (saddles, Mi and M2) in 
the Wp versus Vz plot. The lines correspond to solutions of the variational equations. 
The thresholds for saddle-center bifurcation are at Vz = 0.4278 and 12.2. (b) — (h) 
The equal-energy contours in the Wp versus Wz phase plot for different Vz- A center 
(open circle) S corresponds to a minimum of energy surrounded by closed loops in 
these contour plots. A saddle (solid square). Mi or M2, corresponds to intersection 
of two equal-energy lines. In (d) — (g) the appearance of a saddle and a center for a 
small change of Vz near the thresholds at Vz = 0.4278 and 12.2 is illustrated. In all 
cases ^ = 50 and gdd = —15. 



To perform a linear stability analysis of the axially-symmetric states, we rewrite 
g and (|5|) as US] 

(7) 

+ ^=^^ 2a-adde('^] , (8) 



Xi 
Xs 

±4 



^3, 
1 

0^3 

Jb 1 



2w X\X2 



0^3 

X2 



2a - adde 

1 2N 

2n{x\x\ 



Xx\ 

X2J 



Ciddh 



OCi\ 

X2J i 



(9) 



exp(xi) 

where (xi, X2, X3, X4) = {wp^Wz^Wp^w^)- These equations for widths can be written 
as X = f(x), where x G {xi, X2, X3, X4}. If x^^^ denote the fixed points with x = 0, 
so that f(x*^°^) 0, then the linearization matrix is J = 9f (x)/9x|x^x(o) [16j. An 
examination of eigenvalues of J reveals the nature of stability of the states. The 
eigenvalues come in pairs ±A and lead to exponential growth unless all of them are 



Two-dimensional dipolar Bose-Einstein condensate bright and vortex solitons 



7 



imaginary corresponding to a spectrally stable equilibrium, which is of interest in the 
present context. For the parameters of figure [l] (a), there exist the saddle Mi at 
= (0.9379,0.2147,0,0), and the center S at x^^) = (2.9497,0.5976,0,0), and the 
saddle M2 at x(°) = (22.8943,2.3485,0,0) of which S and M2 are shown in figure § 
(h). The eigenvalues are (±4.8644, ±19.924H) for Mi, (±4.4835z, ±0.2274i) for S, 
and (±0.4751, ±0.0018i) for M2. The center with pairs of pure imaginary eigenvalues 
confirm its stability. 

Using variational equations, we analyze the appearance of axially-symmetric bright 
solitons using the phase plots of \gdd\ versus g for 14 = 2 in figure [3] (a). For \gdd\ in a 
window of critical values, stable bright solitons can be formed. For smaller \gdd\i there is 
too much repulsion and the system expands to infinity and for larger \gdd\^ there is too 
much attraction leading to collapse allowing only unstable stationary states. In figure [3] 
(b) we plot the variational and numerical rms sizes and energies for \gdd\ corresponding 
to the numerical points (^) in figure [3] (a). 

Next, we consider a stable soliton array by mounting tiny DBEC bright solitons 
along the supporting OL sites. Such solitons will be interacting due to long-range 
dipolar interaction. We prepare a stable array of axially-symmetric solitons each 
with g = —gdd = 5 by putting them on alternate sites of OL Vq^ = — 2cos(2z) at 
X = y = 0^ z = ±(2n + l)7r, n = 0, 1, 2, 3. Such an array of bright solitons with small g 
and \gdd\ are stable, whereas those of large g and \gdd\i e.g. the one of figure [l] (a), are 
unstable. In figure [i] (a) and (b) we show the initial array and the the final profile after 
real-time propagation at t = 200. A similar array of anisotropic solitons has a finite life 
and is destroyed at large times. 

Finally, we investigate the collision of two axially-symmetric bright solitons with 




Figure 3. (Color online) (a) The phase plot of \gdd\ versus g from variational analysis 
of the axially-symmetric solitons showing the region of stable solitons for Vz = 2. 
The ^'s denote the numerical points showing the stable-unstable boundary, (b) The 
numerical (n) and variational (v) rms sizes and energy versus g for \gdd\ corresponding 
to the ^'s in (a). 
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Figure 4. (Color online) 3D contour of the stable array of eight axially- symmetric 
solitons each with g = —gdd = 5 on OL Vq£ = — 2cos(2z) at (a) t = and at (b) 
t = 200. The density at contour is 0.001. (c) Contour plot of density 0, t)p of 

two colliding solitons each with g = 5 and gdd = —3 and velocity 1, before, during, 
and after collision. 

g = 5 and gdd = —3 each, placed at x = ±12.8, y = z = Ositt = 0. Each soliton is 
given a velocity of ^ = 1 towards center x = by a phase factor exp(±ix), respectively, 
in the initial wave functions. The collision dynamics is illustrated in figure [i] (c) where 
we show the snapshots of contour plots of density 0, z, t)p at different times. The 
solitons come towards each other, interact at x = and t ^ 12.8 and come out without 
deformation showing their robustness. In this simulation, not only the parameters of 
two individual solitons should lead to a stable state, the combined nonlinearities 2g and 
2gdd should also correspond to a stable configuration in figure |3] (a) to avoid collapse 
during collision. It was demonstrated in \L3\ that, under harmonic confinement, after 
collision at very low velocities, two quasi-2D dipolar BEC solitons may merge together 
to form a single soliton molecule. However, in the present case the solitons appear in a 
narrow window of nonlinearities g and gdd^ as can be seen from figure [3| If two equal 
2D solitons, as in figure [3] (c), coalesce at low velocities, the nonlinearities of the merged 
soliton molecule will be outside the domain of stability in figure [3} Hence, the formation 
of soliton molecule is mostly not possible in the present case. At large velocities the 
solitons of |13j undergo quasi-elastic collision quite similar to the present collision shown 
in figure [3] (c). Also, both quasi-elastic collision at large velocities and merging at low 
velocities of two quasi- ID solitons under transverse harmonic confinement was illustrated 
in [ig. 

To summarize, we studied different types of 2D bright solitons in a DBEC with 
repulsive atomic interaction using the solution of the 3D GP equation. Anisotropic 
stable 2D bright solitons in DBEC are possible on a weak ID OL perpendicular to the 
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polarization direction. Axially-symmetric stable 2D bright and vortex solitons in DBEC 
can be generated on a weak ID OL along the polarization direction when the dipolar 
interaction is tuned to negative values [lOj. In this sign-changed dipolar-interaction 
configuration, bright and vortex solitons are stable due to the long-range attractive 
dipolar interaction in the quasi-2D shape. In the axially-symmetric case, an ID stable 
array of tiny solitons, placed on alternate OL sites, can be made. Such ID array with 
empty OL sites between solitons is of interest in condensed matter physics 
and bears some similarity with stable checkerboard pattern of DBEC on 2D OL with 
empty sites in between [27], both arising due to dipolar interaction. The elastic nature of 
collision of two axially-symmetric solitons is also demonstrated. With present technology 
these stable 2D solitons and their ID arrays can be created and studied in laboratory. 
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